Figure 03 Bile Acids
knitr::opts_chunk$set(
warning = FALSE,
error = TRUE,
echo = TRUE,
fig.width = 10,
fig.height = 7)library("magrittr")
library("data.table")
library("ggplot2"); packageVersion("ggplot2")## [1] '3.3.5'library("patchwork"); packageVersion("patchwork")## [1] '1.1.1'library("ggbeeswarm"); packageVersion("ggbeeswarm")## [1] '0.6.0'theme_set(
theme_bw() +
theme(
panel.grid = element_blank(),
axis.ticks = element_line(size = 0.25),
strip.background = element_blank(),
strip.text.y = element_text(angle = 0)
)
)
scaleColorManualTreatment <-
scale_color_manual(
values =
c(
placebo = "gray50",
wbf10 = "darkblue",
wbf11 = "darkgreen"
)
)
scaleFillManualTreatment <-
scale_fill_manual(
values =
c(
placebo = "gray50",
wbf10 = "darkblue",
wbf11 = "darkgreen"
)
)
studyArmPretty <-
c(placebo = "Placebo",
wbf10 = "WBF-010",
wbf11 = "WBF-011")
scaleColorTreatmentManPretty <-
scale_color_manual(
values =
c(
"Placebo" = "gray50",
"WBF-010" = "darkblue",
"WBF-011" = "darkgreen"
)
)
scaleFillTreatmentManPretty <-
scale_fill_manual(
values =
c(
"Placebo" = "gray50",
"WBF-010" = "darkblue",
"WBF-011" = "darkgreen"
)
)
showBileAcids <- c(
"Cholic acid" = "CA",
"Chenodeoxycholic acid" = "CDCA",
"Ursodeoxycholic acid" = "UDCA")Load data
# Targeted plasma bile acids
tabPlasmaMsOmicsBileAcidsLong <- readRDS(params$tabPlasmaMsOmicsBileAcidsLong)
# "negative" concentration should be set to LoD to avoid logic error
tabPlasmaMsOmicsBileAcidsLong[(Concentration < 0.0), Concentration := LOD]
# Untargeted metabolites change stats (to show the plasma bile acids)
tabMetabolonUntgtWilcox <- readRDS(params$tabMetabolonUntgtWilcox)Fig. 3a Untargeted Bile Acids
Plasma untargeted bile acids summary
# This is both the inclusion set and order it appears in the plots
showConjuFams <- c("UDCA", "CDCA", "CA", "DCA", "LCA")
# Force reset
pBileAcidUntgt <- tabBileAcidUntgt <- NULL
## Untargeted bile acid table
tabBileAcidUntgt <-
tabMetabolonUntgtWilcox %>% copy %>%
# Show the bile acids result in bile-acid figure, omit here to save figure space
.[(subPathway == "Bile Acids")]
# Standardize nomenclature
tabBileAcidUntgt[, Molecule := copy(BIOCHEMICAL)]
tabBileAcidUntgt[, Molecule := gsub("late", "lic acid", Molecule)]
# Define conjugate families
tabBileAcidUntgt[, ConjugateFamily := "Other"]
tabBileAcidUntgt[grep("^([Tt]auro|[Gg]lyco)*[Cc]holic acid", Molecule),
ConjugateFamily := "CA"]
tabBileAcidUntgt[grep("^([Z7]-)*([Kk]eto|[Tt]auro|[Gg]lyco)*[Dd]eoxycholic acid", Molecule),
ConjugateFamily := "DCA"]
tabBileAcidUntgt[grep("^([Tt]auro|[Gg]lyco)*[Cc]hen(o)*deoxycholic acid", Molecule),
ConjugateFamily := "CDCA"]
tabBileAcidUntgt[grep("^([Ii]so)*([Tt]auro|[Gg]lyco)*[Uu]rsodeoxycholic acid", Molecule),
ConjugateFamily := "UDCA"]
tabBileAcidUntgt[grep("[lL]ithocholic acid", Molecule),
ConjugateFamily := "LCA"]
tabBileAcidUntgt[grep("[Hh]yocholic acid", Molecule),
ConjugateFamily := "HCA"]
tabBileAcidUntgt[grep("[Hh]yodeoxycholic acid", Molecule),
ConjugateFamily := "HDCA"]
tabBileAcidUntgt[grep("muri", Molecule, ignore.case = TRUE),
ConjugateFamily := "MCA"]
tabBileAcidUntgt[, OrderBA := "Z"]
tabBileAcidUntgt[(ConjugateFamily %in% c("CA", "CDCA", "HCA")), OrderBA := "Primary"]
tabBileAcidUntgt[(ConjugateFamily %in% c("DCA", "LCA", "UDCA")), OrderBA := "Secondary"]
tabBileAcidUntgt[grep("(keto|iso)", Molecule, ignore.case = TRUE), OrderBA := "Secondary"]
# Upper-case-ify
tabBileAcidUntgt$Molecule <-
tabBileAcidUntgt$Molecule %>%
Hmisc::capitalize()
# Define sort order based on WBF-011
orderUntgtMetabs <-
tabBileAcidUntgt %>%
.[(treatment == "wbf11")] %>%
setorder(-estimateWinGrp) %>%
.$Molecule
tabBileAcidUntgt[, metaboFactor := factor(Molecule, levels = orderUntgtMetabs)]Define base plot
pBileAcidUntgt <-
tabBileAcidUntgt %>% copy %>%
.[(ConjugateFamily %in% showConjuFams)] %>%
.[, conjuFam := factor(ConjugateFamily, levels = showConjuFams)] %>%
# .[(pBwGrp < 0.05)] %>%
ggplot(aes(metaboFactor, estimateWinGrp, color = treatment, fill = treatment)) +
geom_hline(yintercept = 0.0, size = 0.1) +
facet_wrap(
facets = ~conjuFam, nrow = 1,
scales = "free_x", shrink = TRUE, drop = TRUE,
strip.position = "top") +
# Adding this has the effect of centering each panel
geom_blank(mapping = aes(ymin = -confIntLo, ymax = -confIntHi)) +
# The (within-group) log2-ratio + C.I.s
geom_pointrange(
alpha = 1.0,
size = 0.5,
fatten = 3,
shape = 23,
stroke = 0.1,
position = position_dodge(width = 0.6),
mapping = aes(ymin = confIntLo, ymax = confIntHi)
) +
ylab(
expression(paste(
log[2](over(Endpoint,Baseline))
))) +
# xlab("")
scaleColorManualTreatment +
scaleFillManualTreatment +
ggtitle("Plasma, untargeted") +
theme(
legend.position = "none",
panel.grid = element_blank(),
axis.ticks = element_blank(),
axis.ticks.x = element_line(size = 0.25),
# strip.text = element_blank(),
strip.background = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(
size = 10,
angle = 30,
hjust = 1,
vjust = 1))
pBileAcidUntgt
## Show UDCA precursors statistics in WBF-011
tabBileAcidUntgt %>%
.[grep("^(Glyco)*[Ll]ithocholic acid$", Molecule)] %>%
rbind(.,
(tabBileAcidUntgt %>%
.[grep("^(Glyco)*[Cc]henodeoxycholic acid$", Molecule)])
) %>%
.[, .(Molecule, treatment, pBwGrp, estimateWinGrp, pWinGrp)] %>%
setorder(-pWinGrp) %>%
knitr::kable(digits = 2)| Molecule | treatment | pBwGrp | estimateWinGrp | pWinGrp |
|---|---|---|---|---|
| Glycolithocholic acid | wbf10 | 0.25 | 0.02 | 0.98 |
| Chenodeoxycholic acid | wbf10 | 0.50 | -0.02 | 0.89 |
| Glycochenodeoxycholic acid | wbf10 | 0.19 | 0.08 | 0.77 |
| Chenodeoxycholic acid | wbf11 | 0.50 | 0.04 | 0.67 |
| Glycochenodeoxycholic acid | placebo | 0.19 | 0.15 | 0.67 |
| Chenodeoxycholic acid | placebo | 0.50 | -0.12 | 0.58 |
| Glycolithocholic acid | wbf11 | 0.25 | 0.25 | 0.42 |
| Glycolithocholic acid | placebo | 0.25 | -0.39 | 0.39 |
| Glycochenodeoxycholic acid | wbf11 | 0.19 | 0.77 | 0.01 |
## Simplified, untargeted plasma bile acid changes
## Trim down to the less-exotic
## (and usually higher concentration) molecules
filterStringBileAcids <- "(tauro|keto|iso|sulfate|glucuronide)"
pBileAcidUntgt$data <-
pBileAcidUntgt$data[grep(filterStringBileAcids, Molecule,
ignore.case = TRUE, invert = TRUE)]
pBileAcidUntgt
## Table
## Summarize key differences in a table.
pBileAcidUntgt$data[(pBwGrp < 0.05)][(pWinGrp < 0.1)][(treatment == "wbf11")]## treatment BIOCHEMICAL pWinGrp statistic estimateWinGrp
## 1: wbf11 ursodeoxycholate 0.08942025 125 0.5807720
## 2: wbf11 glycoursodeoxycholate 0.04455948 145 0.5085894
## confIntLo confIntHi subPathway estimateBwGrp pBwGrp
## 1: -0.11166223 1.322790 Bile Acids 1.325801 0.016203000
## 2: 0.01318107 1.016929 Bile Acids 1.027953 0.006173199
## Molecule ConjugateFamily OrderBA
## 1: Ursodeoxycholic acid UDCA Secondary
## 2: Glycoursodeoxycholic acid UDCA Secondary
## metaboFactor conjuFam
## 1: Ursodeoxycholic acid UDCA
## 2: Glycoursodeoxycholic acid UDCA# Define withing group marking symbol
pBileAcidUntgt$data[, symbolWinGrp := ""]
pBileAcidUntgt$data[(pBwGrp < 0.05 & pWinGrp < 0.1), symbolWinGrp := "."]
pBileAcidUntgt$data[(pBwGrp < 0.05 & pWinGrp < 0.05), symbolWinGrp := "*"]
# Define between-group (bracket) symbol
pBileAcidUntgt$data[, symbolBwGrp := ""]
pBileAcidUntgt$data[(pBwGrp < 0.05 & pWinGrp < 0.1 & treatment == "wbf11"), symbolBwGrp := "*"]
pBileAcidUntgt$data %>%
.[(pBwGrp < 0.05 & pWinGrp < 0.1)] %>%
.[, .(Molecule, ConjugateFamily, estimateWinGrp, pWinGrp, pBwGrp, symbolWinGrp)] %>%
knitr::kable(digits = 3)| Molecule | ConjugateFamily | estimateWinGrp | pWinGrp | pBwGrp | symbolWinGrp |
|---|---|---|---|---|---|
| Ursodeoxycholic acid | UDCA | 0.581 | 0.089 | 0.016 | . |
| Glycoursodeoxycholic acid | UDCA | 0.509 | 0.045 | 0.006 | * |
## Add brackets
tabBileAcidUntgtBrkt <- NULL
tabBileAcidUntgtBrkt <-
pBileAcidUntgt$data %>% copy %>%
.[(pBwGrp < 0.05 & pWinGrp < 0.1)] %>%
.[, xpos := .I]
pBileAcidUntgt <-
pBileAcidUntgt +
# Bw-grp bracket
ggpubr::geom_bracket(
mapping = aes(
y.position = confIntHi,
label = symbolBwGrp,
xmin = xpos - 0.2,
xmax = xpos + 0.2),
tip.length = c(0.2, 0.02),
bracket.nudge.y = 0.4,
color = "black",
vjust = 0.3,
hjust = 0.5,
data = tabBileAcidUntgtBrkt,
label.size = 4,
) +
# Within-group nominal stat sig
geom_text(
data = tabBileAcidUntgtBrkt,
nudge_y = 0.1,
mapping = aes(label = symbolWinGrp, y = confIntHi, x = xpos + 0.2)
)
pBileAcidUntgt
Fig. 3b Targeted Bile Acids
Event-Wide Table
Compute subject-wise Week12, Baseline paired metrics
loopAndIndicatorVars <-
c("treatment", "quantCategory", "detectCategory",
"OrderBA", "ConjugateFamily", "Molecule")
formulaPivotEventWide <-
loopAndIndicatorVars %>%
c("Subject", .) %>%
paste(collapse = " + ") %>%
paste0(" ~ Event", collapse = "") %>%
as.formula()
tabPlasmaTgtBileAcidsEventWide <- NULL
tabPlasmaTgtBileAcidsEventWide <-
tabPlasmaMsOmicsBileAcidsLong %>%
# pivot wide to ensure consistent handling of subject-event-sample and missing
dcast.data.table(
formula = formulaPivotEventWide,
# formula = Treatment + Subject +
# detectCategory + quantCategory +
# ConjugateFamily + Molecule ~ Event,
fill = NA,
fun.aggregate = mean,
value.var = "Concentration")
# How many different quantitative-accuracy bile acids were measured?
tabPlasmaTgtBileAcidsEventWide %>%
.[(quantCategory == "absolut")] %>%
.$Molecule %>%
unique() %>%
sort()## [1] "12-Ketolithocholic acid" "Chenodeoxycholic acid"
## [3] "Cholic acid" "Deoxycholic acid"
## [5] "Glycochenodeoxycholic acid" "Glycocholic acid"
## [7] "Glycodeoxycholic acid" "Glycoursodeoxycholic acid"
## [9] "Hyodeoxycholic acid" "Isolithocholic acid"
## [11] "Lithocholic acid" "Taurochenodeoxycholic acid"
## [13] "Taurocholic acid" "Taurodeoxycholic acid"
## [15] "Taurolithocholic acid" "Taurolithocholic acid sulfate"
## [17] "Tauroursodeoxycholic acid" "Ursodeoxycholic acid"# Define concise table of Dummy values, set to LoD
tabLoDMsOmics <-
tabPlasmaMsOmicsBileAcidsLong %>%
.[, .(Molecule, LOD)] %>%
unique()
# Add LoD back to each entry in 'wide' table
tabPlasmaTgtBileAcidsEventWide <-
tabLoDMsOmics[tabPlasmaTgtBileAcidsEventWide, on = "Molecule"]
tabPlasmaTgtBileAcidsEventWide[(Baseline < LOD), Baseline := LOD]
tabPlasmaTgtBileAcidsEventWide[(Week12 < LOD), Week12 := LOD]
# Compute change statistics
tabPlasmaTgtBileAcidsEventWide <-
tabPlasmaTgtBileAcidsEventWide %>%
.[, Week12MinusBaseline := Week12 - Baseline]
## Test, "Delta": Week12 - Baseline
tabPlasmaTgtBileAcidsEventWide %>%
.[is.infinite(Week12MinusBaseline)] %>%
nrow()## [1] 0tabPlasmaTgtBaWinGrp <- NULL
tabPlasmaTgtBaWinGrp <-
tabPlasmaTgtBileAcidsEventWide %>%
.[, .(
wilcoxOut = list(
wilcox.test(
x = Week12MinusBaseline,
alternative = "two.sided",
conf.int = TRUE,
mu = 0.0,
paired = FALSE)
)
),
by = loopAndIndicatorVars]
# Extract key results from test
extract_wi = function(wilcoxOut){data.table(
pvalue = wilcoxOut$p.value,
statistic = wilcoxOut$statistic,
estimate = wilcoxOut$estimate,
confIntLo = wilcoxOut$conf.int[1],
confIntHi = wilcoxOut$conf.int[2]
)}
tabPlasmaTgtBaWinGrp <-
tabPlasmaTgtBaWinGrp[, extract_wi(wilcoxOut[[1]]),
by = loopAndIndicatorVars]
setorder(tabPlasmaTgtBaWinGrp, -pvalue, treatment)
## Between-Group: WBF-011 v. Placebo
# Compute wilcoxon
tabPlasmaTgtBaBwGrp <- NULL
# Compute the between-group comparison WBF-011 v. Placebo
tabPlasmaTgtBaBwGrp <-
tabPlasmaTgtBileAcidsEventWide %>%
# .[(numNonMissing > 0)] %>%
# .[(numNonzeroLogRatio > thresholdMinNonzeroLogRatios)] %>%
.[(treatment %in% c("wbf11", "placebo"))] %>%
.[, .(
wilcoxOut = list(
try(expr = {
wilcox.test(
x = Week12MinusBaseline[(treatment == "wbf11")],
y = Week12MinusBaseline[(treatment == "placebo")],
alternative = "greater",
conf.int = TRUE,
paired = FALSE)
}, silent = TRUE)
)
),
by = c(loopAndIndicatorVars[-1])]
# Extract key results from test
extract_wi = function(wilcoxOut){
data.table(
pvalue = wilcoxOut$p.value,
statistic = wilcoxOut$statistic,
estimate = wilcoxOut$estimate,
confIntLo = wilcoxOut$conf.int[1],
confIntHi = wilcoxOut$conf.int[2]
)
}
# Check that test didn't fail and return non-sense (missing or funky data)
tabPlasmaTgtBaBwGrp[, success := inherits(wilcoxOut[[1]], "htest"),
by = c(loopAndIndicatorVars[-1])]
tabPlasmaTgtBaBwGrp <-
tabPlasmaTgtBaBwGrp %>%
.[(success), try({extract_wi(wilcoxOut[[1]])}, silent = TRUE),
by = c(loopAndIndicatorVars[-1])]
setorder(tabPlasmaTgtBaBwGrp, -pvalue)
tabPlasmaTgtBaBwGrp %>% .[(quantCategory == "absolut")] %>% .[(pvalue < 0.25)] %>% tail(10)## quantCategory detectCategory OrderBA ConjugateFamily
## 1: absolut Final results Primary CDCA
## 2: absolut Final results Secondary DCA
## 3: absolut Final results Secondary UDCA
## 4: absolut Final results Secondary DCA
## 5: absolut Final results Secondary UDCA
## 6: absolut Final results Secondary DCA
## Molecule pvalue statistic estimate confIntLo
## 1: Chenodeoxycholic acid 0.19858708 157 0.12638182 -0.106434828
## 2: Taurodeoxycholic acid 0.08146491 172 0.03038607 -0.003110294
## 3: Ursodeoxycholic acid 0.07503215 173 0.01851242 -0.003941349
## 4: Deoxycholic acid 0.06598891 175 0.36320082 -0.011845718
## 5: Glycoursodeoxycholic acid 0.06598891 175 0.04471958 -0.002393418
## 6: Glycodeoxycholic acid 0.05283461 178 0.23311086 -0.001903283
## confIntHi
## 1: Inf
## 2: Inf
## 3: Inf
## 4: Inf
## 5: Inf
## 6: Inf# Join the two comparison tables (within- and between-group)
setnames(tabPlasmaTgtBaWinGrp, "pvalue", "pWinGrp")
setnames(tabPlasmaTgtBaWinGrp, "estimate", "estimateWinGrp")
setnames(tabPlasmaTgtBaBwGrp, "pvalue", "pBwGrp")
setnames(tabPlasmaTgtBaBwGrp, "estimate", "estimateBwGrp")
tabPlasmaTgtWilcox <-
(tabPlasmaTgtBaBwGrp %>%
.[, .(Molecule, estimateBwGrp, pBwGrp)]) %>%
.[tabPlasmaTgtBaWinGrp, on = "Molecule"] %>%
copy()Show UDCA precursors statistics in WBF-011
tabPlasmaTgtWilcox %>%
.[grep("^(Glyco)*[Ll]ithocholic acid$", Molecule)] %>%
.[(quantCategory == "absolut")] %>%
rbind(.,
(tabPlasmaTgtWilcox %>%
.[grep("^(Glyco)*[Cc]henodeoxycholic acid$", Molecule)])
) %>%
.[, .(Molecule, treatment, pBwGrp, estimateWinGrp, pWinGrp)] %>%
setorder(-pWinGrp) %>%
knitr::kable(digits = 2)| Molecule | treatment | pBwGrp | estimateWinGrp | pWinGrp |
|---|---|---|---|---|
| Lithocholic acid | placebo | 0.66 | 0.00 | 0.95 |
| Chenodeoxycholic acid | placebo | 0.20 | -0.02 | 0.90 |
| Glycochenodeoxycholic acid | wbf10 | 0.46 | 0.06 | 0.80 |
| Lithocholic acid | wbf11 | 0.66 | 0.00 | 0.80 |
| Chenodeoxycholic acid | wbf10 | 0.20 | -0.03 | 0.67 |
| Glycochenodeoxycholic acid | placebo | 0.46 | 0.20 | 0.46 |
| Lithocholic acid | wbf10 | 0.66 | 0.02 | 0.44 |
| Chenodeoxycholic acid | wbf11 | 0.20 | 0.09 | 0.29 |
| Glycochenodeoxycholic acid | wbf11 | 0.46 | 0.21 | 0.10 |
Define plotting order, organization.
# Define molecule order based on significance results
showOrder <- showSet <- NULL
showSet <-
tabPlasmaTgtBaBwGrp %>% copy %>%
.[(quantCategory == "absolut")] %>%
# Select the groups to show in main plot
.[(ConjugateFamily %in% showConjuFams)] %>%
# .[(pvalue < 0.5)] %>%
.$Molecule
# Define order of labels (a grouping)
showOrder <-
tabPlasmaTgtWilcox %>% copy %>%
.[showSet, on = "Molecule"] %>%
.[(treatment == "wbf11")] %>%
setorder(-estimateWinGrp) %>%
.$MoleculeTargeted Summary, Delta
# Define table for plot
tabPlasmaTargetedBasePlot <-
tabPlasmaTgtWilcox %>% copy %>%
.[showSet, on = "Molecule"] %>%
# Make order match the untargeted results, via orderUntgtMetabs
.[, bileAcids := factor(Molecule,
levels = c(orderUntgtMetabs, "Lithocholic acid"))] %>%
.[, conjuFam := factor(ConjugateFamily, levels = showConjuFams)]
# Define plot
pPlasmaTargetedBase <-
tabPlasmaTargetedBasePlot %>%
ggplot(aes(x = bileAcids, y = estimateWinGrp)) +
# facet_grid(ConjugateFamily ~ ., scales = "free") +
facet_wrap(~conjuFam, scales = "free", nrow = 1) +
geom_hline(yintercept = 0.0, size = 0.1) +
geom_blank(mapping = aes(ymin = -confIntLo, ymax = -confIntHi)) +
geom_pointrange(
position = position_dodge(width = 0.6),
alpha = 1.0,
size = 0.5,
fatten = 3,
shape = 23,
stroke = 0.1,
mapping = aes(
fill = treatment,
color = treatment,
ymin = confIntLo,
ymax = confIntHi
)
) +
scaleColorManualTreatment +
scaleFillManualTreatment +
# coord_flip() +
ylab("\u0394 Concentration [\u03BCM]") +
ggtitle("Plasma, targeted") +
theme(
legend.position = "none",
panel.grid = element_blank(),
axis.ticks.x = element_line(size = 0.25),
# panel.border = element_blank(),
strip.background = element_blank(),
axis.title.y = element_text(size = 14, angle = 90, vjust = 0.5),
axis.title.x = element_blank(),
axis.text.x = element_text(
size = 10,
angle = 30,
hjust = 1,
vjust = 1))
pPlasmaTargetedBase
## Simplified, targeted plasma bile acid changes
## Trim down to the less-exotic
## (and usually higher concentration) molecules
pPlasmaTargetedBase$data <-
pPlasmaTargetedBase$data[grep(filterStringBileAcids, Molecule,
ignore.case = TRUE, invert = TRUE)]
pPlasmaTargetedBase
## Table
## Summarize key differences in a table.
pPlasmaTargetedBase$data[(pBwGrp < 0.15)][(pWinGrp < 0.3)][(estimateWinGrp > 0.0)]## Molecule estimateBwGrp pBwGrp treatment quantCategory
## 1: Ursodeoxycholic acid 0.01851242 0.07503215 wbf11 absolut
## 2: Glycoursodeoxycholic acid 0.04471958 0.06598891 wbf11 absolut
## detectCategory OrderBA ConjugateFamily pWinGrp statistic estimateWinGrp
## 1: Final results Secondary UDCA 0.2012040 104 0.01232496
## 2: Final results Secondary UDCA 0.2412529 125 0.01250471
## confIntLo confIntHi bileAcids conjuFam
## 1: -0.006103741 0.03491730 Ursodeoxycholic acid UDCA
## 2: -0.007507829 0.03130284 Glycoursodeoxycholic acid UDCA# Define withing group marking symbol
pPlasmaTargetedBase$data[, symbolWinGrp := ""]
pPlasmaTargetedBase$data[(pBwGrp < 0.15 & pWinGrp < 0.3 & estimateWinGrp > 0.0), symbolWinGrp := "."]
# pPlasmaTargetedBase$data[(pBwGrp < 0.05 & pWinGrp < 0.05), symbolWinGrp := "*"]
# Define between-group (bracket) symbol
pPlasmaTargetedBase$data[, symbolBwGrp := ""]
pPlasmaTargetedBase$data[(pBwGrp < 0.15 & pWinGrp < 0.3 &
estimateWinGrp > 0.0 &
treatment == "wbf11"), symbolBwGrp := "."]
pPlasmaTargetedBase$data %>%
.[(pBwGrp < 0.15 & pWinGrp < 0.3 &
estimateWinGrp > 0.0)] %>%
.[, .(Molecule, ConjugateFamily, estimateWinGrp, pWinGrp, pBwGrp, symbolWinGrp)] %>%
knitr::kable(digits = 3)| Molecule | ConjugateFamily | estimateWinGrp | pWinGrp | pBwGrp | symbolWinGrp |
|---|---|---|---|---|---|
| Ursodeoxycholic acid | UDCA | 0.012 | 0.201 | 0.075 | . |
| Glycoursodeoxycholic acid | UDCA | 0.013 | 0.241 | 0.066 | . |
## Add brackets
tabBileAcidTgtBrkt <- NULL
tabBileAcidTgtBrkt <-
pPlasmaTargetedBase$data %>% copy %>%
.[(symbolBwGrp != "")] %>%
.[, xpos := .I]
pPlasmaTargeted <-
pPlasmaTargetedBase +
# Bw-grp bracket
ggpubr::geom_bracket(
mapping = aes(
y.position = confIntHi,
label = symbolBwGrp,
xmin = xpos - 0.2,
xmax = xpos + 0.2),
tip.length = c(0.2, 0.1),
bracket.nudge.y = 0.04,
color = "black",
vjust = 0.15,
hjust = 0.5,
data = tabBileAcidTgtBrkt,
label.size = 6,
) +
# Within-group nominal stat sig
# (nothing: within-group sig weaker than the untargeted log-ratio data)
geom_text(
data = tabBileAcidTgtBrkt[(pWinGrp < 0.1)],
nudge_y = 0.01,
mapping = aes(
color = treatment,
label = symbolWinGrp,
y = confIntHi,
x = xpos + 0.2)
)
pPlasmaTargeted
Fig. 3c Total bile acids
For completeness, sum up total bile acids, show whether it is within healthy range for all subjects.
- Annotated the reference ranges.
- Annotate the medians within group, and grand at each timepoint.
studyArmDodgeWidth = 0.8
thresholdMaxTotalPlot = 15
tabShowBileAcids <- pPlasmaTotalBileAcids <- NULL
tabShowBileAcids <-
tabPlasmaMsOmicsBileAcidsLong %>% copy %>%
.[(
detectCategory == "Final results" &
quantCategory == "absolut")] %>%
.[, .(totalBileAcids = sum(Concentration, na.rm = TRUE)),
by = .(treatment, Subject, Event)] %>%
.[, Treatment := studyArmPretty[treatment]]
tabPlotTotalBileAcids <-
tabShowBileAcids %>% copy %>%
.[(totalBileAcids > thresholdMaxTotalPlot),
totalBileAcids := Inf]
pPlasmaTotalBileAcids <-
tabPlotTotalBileAcids %>%
ggplot(aes(Treatment, totalBileAcids, color = Treatment, fill = Treatment)) +
facet_wrap(~Event, strip.position = "top") +
# Normal range, 0 - 10 uM
# Moderate cholestasis ~10 - 40 uM
geom_rect(
data = data.frame(Treatment = "dummy", Subject = "SS"),
color = NA,
fill = "darkred",
alpha = 0.15,
inherit.aes = FALSE,
xmin = -Inf,
xmax = Inf,
ymin = 10,
ymax = Inf
) +
# severe cholestasis
# >= 40 uM
# (not shown)
#
# Baseline grand median
geom_hline(
yintercept = tabPlotTotalBileAcids[(Event == "Baseline")] %>%
.$totalBileAcids %>%
median(),
size = 0.25,
linetype = 1,
alpha = 0.65,
) +
geom_beeswarm(
data = tabPlotTotalBileAcids[is.finite(totalBileAcids)],
shape = 21,
dodge.width = studyArmDodgeWidth,
# mapping = aes(color = Treatment, fill = Treatment),
groupOnX = TRUE,
cex = 3.5,
size = 1,
stroke = 0.1,
alpha = 0.75
) +
# Annotate any points that get off the top of the screen
ggrepel::geom_text_repel(
position = position_dodge(width = studyArmDodgeWidth),
data = tabShowBileAcids[(totalBileAcids > thresholdMaxTotalPlot)],
mapping = aes(
y = Inf,
label = round(totalBileAcids, digits = 1)
),
min.segment.length = 0,
size = 2
) +
# group-median crossbar
stat_summary(
geom = "crossbar",
fatten = 1.3,
width = 0.5,
position = position_dodge(width = studyArmDodgeWidth),
fun.data = function(x){
thisMedian = median(x, na.rm = TRUE)
return(
data.frame(
y = thisMedian,
ymin = thisMedian,
ymax = thisMedian
)
)
}) +
scaleColorTreatmentManPretty +
scaleFillTreatmentManPretty +
guides(
color = guide_legend(title="Study Arm")
) +
ylab("Total bile acids [\u03BCM]") +
theme(
panel.grid = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_line(size = 0.25),
axis.text.x = element_text(angle = 30, size = 8, hjust = 1),
legend.position = "none"
)
pPlasmaTotalBileAcids
Fig. 3e Monoculture \(\Delta\) Bile Acids
Load data
tabStrainCultureMsOmicsBileAcids <-
readRDS(params$tabStrainCultureMsOmicsBileAcids)Check NIC
tabStrainCultureMsOmicsBileAcids %>%
.[(inoculated == "FALSE")] %>%
.[(sampleType == "supernatant")] %>%
.[(quantCategory == "absolut")] %>%
ggplot(aes(timepoint, Concentration, color = ConjugateFamily)) +
facet_grid(cols = vars(moleculeAdded), rows = vars(run)) +
geom_hline(yintercept = 50, linetype = 3, size = 0.25) +
geom_point() +
theme(axis.ticks.y = element_line(size = 0.25)) +
ggtitle("NIC only", "Looks good. No changes. Confirms design.")
Check growth-typical conditions (no added bile acids) controls. (note: this control was not included in CBUT run)
tabStrainCultureMsOmicsBileAcids %>% copy %>%
.[(Condition == "GrowthCtl")] %>%
.[(sampleType == "supernatant")] %>%
.[(quantCategory == "absolut")] %>%
ggplot(aes(timepoint, Concentration, color = ConjugateFamily)) +
facet_grid(cols = vars(Strain)) +
ylim(0, 50) +
geom_path(mapping = aes(group = paste(Condition, Molecule))) +
geom_point() +
theme(
axis.ticks.y = element_line(size = 0.25),
legend.position = "bottom") +
ggtitle("Control: No primary bile acids amended",
"Minimal changes, low concentrations. Confirms design.")
Note: Should include only “absolut” concentration values in the figure-of-merit.
Effective volume, umoles
Compute the effective-volume that combines the cell pellet and supernatant measurements at t_final.
tabMonocultureDeltaSummary <-
tabStrainCultureMsOmicsBileAcids %>% copy %>%
.[(inoculated == "TRUE")] %>%
.[(Condition != "GrowthCtl")] %>%
.[(media == "pyg")] %>%
.[(quantCategory == "absolut")] %>%
# This sums together cell-pellet and supernatant at t_final,
# on both micromoles and effective_volume values
dcast.data.table(
formula = inoculated + moleculeAdded + Strain +
Condition +
ConjugateFamily + OrderBA + Molecule ~ timepoint,
value.var = c("micromoles", "effectiveVolume"),
fun.aggregate = sum,
na.rm = TRUE
) %>%
# Compute effective concentrations at zero, final
.[, conc0 := (1000 * micromoles_T0 / effectiveVolume_T0)] %>%
.[, concf := (1000 * micromoles_Tfinal / effectiveVolume_Tfinal)] %>%
# Deltas
.[, deltaMicroMoles := micromoles_Tfinal - micromoles_T0] %>%
.[, deltaConc := concf - conc0]Show each timepoint separately.
tabPlotBothTimepoints <-
tabStrainCultureMsOmicsBileAcids %>%
.[!is.na(Strain)] %>%
.[(inoculated == "TRUE")] %>%
.[(Condition != "GrowthCtl")] %>%
.[(media == "pyg")] %>%
.[(quantCategory == "absolut")] %>%
# Sum
.[, .(micromoles = sum(micromoles)),
by = .(Strain, moleculeAdded, Condition,
ConjugateFamily, Molecule, timepoint)]
# Only show molecules that have micromoles > threshold
showMolecules <- tabPlotBothTimepoints[(micromoles > 1e-3)]$Molecule %>% unique()
tabPlotBothTimepoints %>%
.[showMolecules, on = "Molecule"] %>%
ggplot(aes(timepoint, micromoles,
# shape = sampleType,
color = ConjugateFamily)) +
facet_grid(
cols = vars(moleculeAdded),
rows = vars(Strain)) +
geom_path(mapping = aes(group = paste(Strain, Condition, Molecule))) +
geom_point() +
ggtitle("Inoculated only", "Replicates similar. Confirms design.")
Summarize as effective concentration deltas.
only abs(\u0394) > 1\u03BCM"
CBUT monoculture main bile acid changes following growth in rich media amended with CA, CDCA, or sterile water (control)
pMonocultureDeltaSummary <- NULL
strainOrder <- c("CBUT", "EHAL", "CBEI", "AMUC", "BINF")
tabDeltaSummary <-
tabMonocultureDeltaSummary %>% copy %>%
.[(Molecule %in% names(showBileAcids))] %>%
# Pretty plotting
.[, Amendment := paste0("+50 \u03BCM ", moleculeAdded)] %>%
# Omit the vanishing trace bile acids from this plot
.[showMolecules, on = "Molecule", nomatch = 0] %>%
# set factor to order the stain facets
.[, strainFac := factor(Strain, levels = strainOrder)] %>%
# limit negative-delta outlier effect on plot limits
.[(deltaConc < -75), deltaConc := -Inf]
pMonocultureDeltaSummary <-
tabDeltaSummary %>%
.[!is.na(Strain)] %>%
ggplot(
mapping = aes(
x = Molecule,
y = deltaConc,
color = OrderBA,
fill = OrderBA)
) +
scale_y_continuous(breaks = c(-50, -25, 0, 25)) +
facet_grid(
rows = vars(Amendment),
cols = vars(strainFac)
) +
geom_hline(yintercept = 0.0, size = 0.25, color = "black") +
# The bars
stat_summary(
width = 0.5,
alpha = 0.7,
geom = "col",
color = NA,
fun.data = function(x){
data.frame(y = median(x))
}) +
# The individual replicate values
geom_point(
size = 1.5,
stroke = 0,
alpha = 0.7,
position = position_beeswarm(cex = 2.5, groupOnX = TRUE)
) +
guides(
colour = guide_legend("Bile Acid Rank"),
fill = guide_legend("Bile Acid Rank")) +
scale_color_brewer(palette = "Set1") +
scale_fill_brewer(palette = "Set1") +
# scale_y_continuous(breaks = c(-40, -20, 0, 20)) +
ylab("\u0394 Concentration [\u03BCM]") +
# ggtitle("\u0394 Concentration [\u03BCM]") +
theme(
plot.title = element_blank(),
strip.text.x = element_text(size = 12),
strip.text.y = element_text(size = 10, hjust = 0),
legend.position = "none",
panel.grid.major.y = element_line(size = 0.1, color = "gray"),
axis.ticks = element_blank(),
axis.title.y = element_text(size = 12),
axis.text.x = element_text(angle = 30, hjust = 1, size = 10),
axis.title.x = element_blank(),
)
# pMonocultureDeltaSummaryFig. 3e Monoculture \(\Delta\) Summary
pMonocultureDeltaSummary## Error in .calculateSwarmUsingC(xy$y, dsize = ysize * cex, gsize = xsize * : NA/NaN/Inf in foreign function call (arg 1)
For plot, make a version that only shows the targeted assay positive confirmation (CBUT). Can leave as “data not shown” for now for the others, since they are all UDCA-negligible.
pMonocultureDeltaSummaryCbutOnly <- copy(pMonocultureDeltaSummary)
pMonocultureDeltaSummaryCbutOnly$data <-
pMonocultureDeltaSummary$data[(Strain == "CBUT")] %>%
copy()
# Flip the facet grid for better vertical space usage
pMonocultureDeltaSummaryCbutOnly <-
pMonocultureDeltaSummaryCbutOnly +
facet_grid(
cols = vars(Amendment),
rows = vars(strainFac)
)
pMonocultureDeltaSummaryCbutOnly
Fig. 3d Strain Culture Pilot, Metabolon
Strain culture rich media + primary bile acids (CA + CDCA) pilot. Metabolon untargeted data.
- 50 uM CA amended
- 50 uM CDCA amended
- No bile acid amendment for AMUC
tabStrainCultureMetabolonPivWideType <-
# Uses Metabolon `OrigScale` rather than the rescaled, imputed values.
# Has benefit of explicit missing values, and consistent peak-area relationship
readRDS(params$tabStrainCultureMetabolonPivWideType)
tabStrainCultureMetabolon <-
readRDS(params$tabStrainCultureMetabolon)
# Subset to bile acids on the 'sampletype wide' table
vecBileAcidIds <-
tabStrainCultureMetabolon[grep("Bile", `SUB PATHWAY`)] %>%
.$COMP_ID %>% unique()
# Define bile acids table
tabStrainCultureMetabolonBileAcids <- NULL
tabStrainCultureMetabolonBileAcids <-
tabStrainCultureMetabolonPivWideType %>% copy %>%
.[vecBileAcidIds, on = "COMP_ID"] %>%
# One is a very small value on this `OrigScale` scale
.[is.na(CellPellet), CellPellet := 1.0] %>%
.[is.na(BlankMedia), BlankMedia := 1.0] %>%
.[is.na(Supernatant), Supernatant := 1.0] %>%
# Compute log-ratios
.[, CpOverBlank := log10(CellPellet / BlankMedia)] %>%
.[, SuOverBlank := log10(Supernatant / BlankMedia)] %>%
# Melt (pivot long) so that can be combined
melt.data.table(
id.vars = c("COMP_ID", "Strain"),
measure.vars = c("CpOverBlank", "SuOverBlank"),
variable.name = "SampleType",
value.name = "log10Ratio") %>%
.[, SampleType := gsub("SuOverBlank", "Supernatant", SampleType)] %>%
.[, SampleType := gsub("CpOverBlank", "CellPellet", SampleType)]
tabStrainCultureMetabolonBileAcids <-
tabStrainCultureMetabolonBileAcids %>%
# Add back annotations
(tabStrainCultureMetabolon %>%
.[, .(COMP_ID, BIOCHEMICAL, `SUPER PATHWAY`, `SUB PATHWAY`,
MASS, CAS, KEGG, `Group HMDB`)] %>%
unique())[., on = "COMP_ID", nomatch = 0]
# Standardize nomenclature, and categorize conjugate family
tabStrainCultureMetabolonBileAcids[, Molecule := copy(BIOCHEMICAL)]
tabStrainCultureMetabolonBileAcids[, Molecule := gsub("late", "lic acid", Molecule)]
tabStrainCultureMetabolonBileAcids[, ConjugateFamily := "Other"]
# Define conjugate families
tabStrainCultureMetabolonBileAcids[grep("^([Tt]auro|[Gg]lyco)*(3-dehydro)*[Cc]holic acid", Molecule),
ConjugateFamily := "CA"]
tabStrainCultureMetabolonBileAcids[grep("^([Z7]-)*([Kk]eto|[Tt]auro|[Gg]lyco)*[Dd]eoxycholic acid", Molecule),
ConjugateFamily := "DCA"]
tabStrainCultureMetabolonBileAcids[grep("^([Tt]auro|[Gg]lyco)*(3-dehydro)*[Cc]hen(o)*deoxycholic acid", Molecule),
ConjugateFamily := "CDCA"]
tabStrainCultureMetabolonBileAcids[grep("^([Ii]so)*([Tt]auro|[Gg]lyco)*[Uu]rsodeoxycholic acid", Molecule),
ConjugateFamily := "UDCA"]
tabStrainCultureMetabolonBileAcids[grep("^([Tt]auro|[Gg]lyco)*[Uu]rsocholic acid", Molecule),
ConjugateFamily := "UCA"]
tabStrainCultureMetabolonBileAcids[grep("[lL]ithocholic acid", Molecule),
ConjugateFamily := "LCA"]
tabStrainCultureMetabolonBileAcids[grep("[Hh]yocholic acid", Molecule),
ConjugateFamily := "HCA"]
tabStrainCultureMetabolonBileAcids[grep("[Hh]yodeoxycholic acid", Molecule),
ConjugateFamily := "HDCA"]
tabStrainCultureMetabolonBileAcids[grep("muri", Molecule, ignore.case = TRUE),
ConjugateFamily := "MCA"]
tabStrainCultureMetabolonBileAcids[, OrderBA := "Z"]
tabStrainCultureMetabolonBileAcids[(ConjugateFamily %in% c("CA", "CDCA", "HCA")),
OrderBA := "Primary"]
tabStrainCultureMetabolonBileAcids[(ConjugateFamily %in%
c("DCA", "LCA", "UDCA", "UCA", "HDCA")),
OrderBA := "Secondary"]
tabStrainCultureMetabolonBileAcids[grep("(keto|iso)", Molecule, ignore.case = TRUE),
OrderBA := "Secondary"]
# For interpretive organization, it helps to group these together, last
tabStrainCultureMetabolonBileAcids[grep("3-dehydro", Molecule),
ConjugateFamily := "z3dehydro"]
tabStrainCultureMetabolonBileAcids[grep("3-dehydro", Molecule),
OrderBA := "Secondary"]
# Upper-case-ify
tabStrainCultureMetabolonBileAcids$Molecule <-
tabStrainCultureMetabolonBileAcids$Molecule %>%
Hmisc::capitalize()
# Define order of metabolites
orderStrainCultureMetabolonBileAcids <-
tabStrainCultureMetabolonBileAcids %>%
.[, .(OrderBA, ConjugateFamily, Molecule)] %>%
unique() %>%
setorder(-OrderBA, -ConjugateFamily, -Molecule) %>%
# show()
.$Molecule
# Define strain culture pilot summary heatmap
pStrainCultureMetabolonBileAcids <-
tabStrainCultureMetabolonBileAcids %>% copy %>%
.[, MoleculeFac := factor(Molecule,
levels = orderStrainCultureMetabolonBileAcids)] %>%
.[, SampleTypePretty := c(CellPellet = "Cell Pellet",
Supernatant = "Supernatant")[SampleType]] %>%
ggplot(aes(SampleTypePretty, MoleculeFac, fill = log10Ratio)) +
# facet_grid(ConjugateFamily ~ Strain, shrink = FALSE, drop = TRUE, scales = "free") +
facet_grid( ~ Strain, shrink = FALSE, drop = TRUE, scales = "free") +
geom_raster() +
scale_fill_gradient2() +
# Adjust the color mapping legend
guides(
fill = guide_colourbar(
title = expression(paste(
log[10](over(Specimen,Medium))
)),
# nbin = 5,
barwidth = 13,
barheight = 0.5)
) +
theme(
legend.position = "bottom",
legend.justification = "left",
legend.title = element_text(size = 8, vjust = 0.2),
axis.title = element_blank(),
axis.text.x = element_text(angle = 30, size = 10, hjust = 1)
)
pStrainCultureMetabolonBileAcids
pStrainCultureMetabolonBileAcids$data <-
# Drop AMUC, because no primary bile acids amended
pStrainCultureMetabolonBileAcids$data %>%
.[(Strain != "AMUC")]
pStrainCultureMetabolonBileAcids
Build Figure 3
fig2Layout <- "
AAAAAAAA#
BBBBBBBCC
DDDD##EEE
"
sizeAxisTitles <- 7
sizeAxisText <- 6
pFig3 <- NULL
pFig3 <-
# Untargeted, Log-Ratio
(
pBileAcidUntgt +
theme(
panel.spacing = unit(8, "pt"),
axis.ticks.y = element_line(size = 0.25),
plot.title = element_text(size = sizeAxisTitles),
strip.text = element_text(size = sizeAxisTitles),
axis.title.x = element_blank(),
axis.text.x = element_text(size = sizeAxisText),
axis.text.y = element_text(size = sizeAxisText),
axis.title.y = element_text(
size = sizeAxisTitles, angle = 0, vjust = 0.5,
margin = margin(r = -100, unit = "pt"))
)
) +
# Targeted, Delta
(
pPlasmaTargeted +
theme(
panel.spacing = unit(1, "pt"),
plot.title = element_text(size = sizeAxisTitles),
strip.text = element_text(size = sizeAxisTitles),
axis.title.x = element_blank(),
axis.text.x = element_text(size = sizeAxisText),
axis.text.y = element_text(size = sizeAxisText),
axis.title.y = element_text(
size = sizeAxisTitles, angle = 90, vjust = 0.5,
margin = margin(r = -100, unit = "pt"))
)
) +
# Targeted, total bile acids
(pPlasmaTotalBileAcids +
theme(
panel.spacing = unit(1, "pt"),
plot.margin = margin(0, unit='pt'),
plot.title = element_text(size = sizeAxisTitles),
strip.text = element_text(size = sizeAxisTitles),
axis.title.x = element_blank(),
axis.title.y = element_text(size = sizeAxisTitles),
axis.text.x = element_text(size = sizeAxisText),
axis.text.y = element_text(size = sizeAxisText)
)
) +
# In vitro monoculture pilot summary heatmap
(
pStrainCultureMetabolonBileAcids +
# Adjust the color mapping legend
guides(
fill = guide_colourbar(
title = expression(paste(
log[10](over(Specimen,Medium))
)),
barwidth = 10,
barheight = 0.5)
) +
theme(
panel.spacing = unit(1, "pt"),
plot.title = element_text(size = sizeAxisTitles),
strip.text = element_text(size = sizeAxisTitles),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(size = sizeAxisText),
axis.text.y = element_text(size = sizeAxisText),
legend.text = element_text(size = sizeAxisText),
legend.title = element_text(size = sizeAxisText, vjust = 0.2),
legend.margin = margin(l = -80, r = 200, t = -25, unit='pt')
)
) +
(pMonocultureDeltaSummaryCbutOnly +
theme(
plot.margin = margin(0, unit='pt'),
axis.ticks.y = element_line(size = 0.25),
axis.ticks.x = element_line(size = 0.25),
plot.title = element_blank(),
strip.text.x = element_text(size = sizeAxisTitles),
strip.text.y = element_text(size = sizeAxisTitles),
axis.title.x = element_blank(),
axis.title.y = element_text(size = sizeAxisTitles),
axis.text.x = element_text(size = sizeAxisText),
axis.text.y = element_text(size = sizeAxisText)
)
) +
plot_annotation(tag_levels = 'a')
# render prototype in report
pFig3 +
plot_layout(design = fig2Layout) &
theme(plot.tag = element_text(size = 8, face = "bold"))
figHeight = 180
ggsave(
"Figure-03.png",
pFig3 +
plot_layout(design = fig2Layout) &
theme(plot.tag = element_text(size = 8, face = "bold")),
width = 10.5, height = 10)
ggsave("Figure-03.pdf",
pFig3 +
plot_layout(design = fig2Layout) &
theme(plot.tag = element_text(size = 8, face = "bold", family = "Sans")),
device = cairo_pdf,
width = 170,
dpi = 300,
height = figHeight,
units = "mm")